
******************************************************************************************************************************
**********************************Figure S.4**********************************************************************************
******************************************************************************************************************************

*************************************************************bound test*******

use point,clear
*sij1
forvalues i=1/14{
local falpha`i'=alpha`i'
}
local fbeta=0.5308357
keep dgp
append using flexlogit_uij
drop if ID==.
drop dgp

gen xbeta1 = `falpha1'*x1+`falpha2'*x2+`falpha3'*x3+`falpha4'*x4+`falpha5'*x5+`falpha6'*x6+`fbeta'*z
forvalues i=7/14 {
replace xbeta1=xbeta1+`falpha`i''*x`i'
}

gen expuij=exp(uij)
bys ID: egen sumexpuij = sum(expuij)
gen sij=expuij/sumexpuij

*Now we know the full utility function

*Let's loop through and do R simulations
*In each simulation, compute U and VU for every good
gen othercount = 0
gen bestbothcount = 0

local R = 1000
forvalues r = 1/`R' {
display "Draw is `r'"
set seed 400`r'
qui gen gumbel = -ln(-ln(runiform())) 
qui gen VU = xbeta1+gumbel-`fbeta'*z
qui gen U = xbeta1+gumbel
*For upper bound, compute:
*P(U_ik > U_ij and VU_ik > VU_ij for some k)
*i.e. does some other good have higher utility and higher visible utility?
*Sort by utility, check each good with higher utility and see if it has higher visible utility
gsort ID -U
by ID: gen Urank = _n
qui gen otherbetter = 0
qui replace otherbetter = 1 if Urank == 2 & VU[_n-1] > VU
qui replace otherbetter = 1 if Urank == 3 & ((VU[_n-1]>VU) | (VU[_n-2]>VU))
qui replace otherbetter = 1 if Urank == 4 & ((VU[_n-1]>VU) | (VU[_n-2]>VU)| (VU[_n-3]>VU))
qui replace otherbetter = 1 if Urank == 5 & ((VU[_n-1]>VU) | (VU[_n-2]>VU)| (VU[_n-3]>VU)| (VU[_n-4]>VU))
qui replace otherbetter = 1 if Urank == 6 & ((VU[_n-1]>VU) | (VU[_n-2]>VU)| (VU[_n-3]>VU)| (VU[_n-4]>VU)| (VU[_n-5]>VU))
qui replace otherbetter = 1 if Urank == 7 & ((VU[_n-1]>VU) | (VU[_n-2]>VU)| (VU[_n-3]>VU)| (VU[_n-4]>VU)| (VU[_n-5]>VU)| (VU[_n-6]>VU))
qui replace otherbetter = 1 if Urank == 8 & ((VU[_n-1]>VU) | (VU[_n-2]>VU)| (VU[_n-3]>VU)| (VU[_n-4]>VU)| (VU[_n-5]>VU)| (VU[_n-6]>VU)| (VU[_n-7]>VU))
qui replace otherbetter = 1 if Urank == 9 & ((VU[_n-1]>VU) | (VU[_n-2]>VU)| (VU[_n-3]>VU)| (VU[_n-4]>VU)| (VU[_n-5]>VU)| (VU[_n-6]>VU)| (VU[_n-7]>VU)| (VU[_n-8]>VU))
qui replace otherbetter = 1 if Urank == 10 & ((VU[_n-1]>VU) | (VU[_n-2]>VU)| (VU[_n-3]>VU)| (VU[_n-4]>VU)| (VU[_n-5]>VU)| (VU[_n-6]>VU)| (VU[_n-7]>VU)| (VU[_n-8]>VU)| (VU[_n-9]>VU))
qui replace othercount = othercount + otherbetter
*For lower bound, compute:
*P(U_ij > U_ik and VU_ij > VU_ik for all k)
*i.e. is good j highest in terms of both utility and visible utility?
gsort ID -VU
by ID: gen VUrank = _n
gen bestboth = (Urank == 1)*(VUrank == 1)
replace bestbothcount = bestbothcount + bestboth
drop Urank VUrank bestboth otherbetter VU U gumbel
}


gen Pother = othercount / `R'
gen Pbestboth = bestbothcount / `R'


gen UB = 1-Pother
gen LB = Pbestboth
sort LB
gen xaxis = _n

*output data and graph in matlab
save boundtest,replace

xtile xquant = xaxis, nq(100)
collapse (mean) UB LB sij, by(xquant)
rename (xquant UB LB) (lbquant ub lb)
outsheet lb ub sij lbquant using graphdata.csv, comma replace





